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INTRODUCTION 


During  the  past  decade  and  a  half,  a  great  deal  of  interest  has  arisen  in  the  development 
of  reliable,  robust,  and  efficient  software  for  the  numerical  solution  of  partial  differential 
equations.  The  origin  of  this  interest  is  obvious.  Partial  differential  equations  are  used  to 
mathematically  model  physical  reality  for  many  scientific  disciplines.  ITie  complexity  of  these 
equations  usually  require  that  numerical  methods  be  employed  in  order  to  obtain  solutions. 

Most  scientists  and  engineers  prefer  to  use  existing  software  for  this  process  in  order  to  avoid 
developing  their  own  numerical  techniques.  In  this  way,  they  can  concentrate  on  the 
interpretation  of  the  computed  results. 

Before  this  time,  even  though  the  interest  in  solving  partial  differential  equations 
numerically  existed,  computers  were  not  up  to  the  task  of  supporting  general  purpose  packages 
along  the  lines  of  what  existed  for  ordinary  differential  equations  (cf.,  eg.,  Gear  (ref  1)  or  Petzold 
(ref  2)).  Due  to  restricted  computing  power,  the  tendency  was  to  write  problem  specific 
computer  codes.  However,  advances  in  computer  hardware,  computer  systems,  and  numerical 
algorithms  led  to  the  feasibility  of  developing  general  purpose  partial  differential  equation  solvers 
(cf.,  e.g.,  Babuska,  Chandra,  and  Flaherty  (ref  3)  and  Flaherty,  Paslow,  Shephard,  and 
Vasilakis  (ref  4)  for  a  sampling). 

The  major  contribution  that  the  field  of  numerical  analysis  made  to  this  endeavor  was  the 
development  of  adaptive  numerical  methods  (refs  3,4).  An  adaptive  numerical  method  is  a 
numerical  method  which  automatically  adjusts  the  discretization  of  the  problem  domain  and/or 
the  order  of  the  method  while  the  solution  of  the  problem  is  being  determined.  This  adjustment 
is  performed  so  as  to  optimize  the  solution  process  in  some  way.  This  optimization  usually  takes 
the  form  of  obtaining  a  solution  with  as  coarse  a  discretization  level  and/or  as  low  an  order  as 
possible  for  a  given  error  tolerance.  This  criterion  is  not  applied  globally,  however,  but  in  a  local 
manner.  The  result  is  that  the  same  discretization  and/or  order  is  not  used  uniformly  throughout 
the  solution  process,  as  is  the  case  with  ordinary  numerical  methods.  Rather,  finer 
discretizations  and/or  higher  order  methods  are  used  in  regions  of  the  problem  domain  where 
numerical  approximations  are  more  difficult  to  make,  and  coarser  discretizations  and/or  lower 
orders  are  used  elsewhere.  In  this  way,  a  reliable  numerical  solution  is  found  as  efficiently  as 
possible,  in  terms  of  computer  resources,  and  with  minimal  user  interface. 

Even  with  the  advances  in  computer  science,  electrical  engineering,  and  numerical 
analysis,  it  is  still  not  practical  to  write  a  software  package  that  attempts  to  solve  all  possible 
types  of  partial  differential  equations.  At  the  very  least,  the  three  basic  types  of  partial 
ffifferential  equations:  elliptic,  hyperbolic,  and  parabolic,  demand  three  different  numerical 
approaches  in  order  to  compute  reliable  solutions  (cf.,  Lapidus  and  Finder  (ref  5)).  The  typical 
general  purpose  partial  differential  equation  solver  is  still  going  to  be  limited  in  scope  in  some 
way. 


This  is  the  first  in  a  series  of  four  reports  whose  overall  purpose  is  to  describe  an 
adaptive  finite  element  method  (AFEM)  for  solving  systems  of  parabolic  partial  differential 
equations  of  the  form 
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«,  +  f{x,t,u,u^  =  [DiXyt,u)uX  ,a<x<b,t>0. 


(1) 


subject  to  prescribed  initial  conditions  and  linear  separated  boundary  conditions.  The  scalar 
variables  x  and  t  represent  spatial  and  temporal  coordinates  and  denote  partial  differentiation 
when  used  as  subscripts;  a  and  b  are  scalar  constants;  u(x,t)  and  f(x,t,u,uj  are  M-vectors;  and 
D(x,t,u)  is  an  MxM  matrix.  The  ultimate  goal  of  AFEM  is  to  determine  an  approximate  solution 
to  Eq.  (1)  to  within  a  user  prescribed  error  tolerance. 

Equations  of  the  form  of  Eq.  (1)  arise  in  many  applications.  When  /  is  identically  zero, 
Eq.  (1)  is  a  system  of  nonlinear  heat  or  diffusion  equations.  They  model  physical  circumstances 
as  diverse  as  the  temperature  in  a  melting  solid  (cf.,  Friedman  (ref  6))  and  bacterial  motion  (cf., 
Keller  and  Odell  (ref  7)).  If/ is  a  function  of  u  only,  then  Eq.  (1)  is  a  system  of  reaction- 
diffusion  equations.  Such  equations  have  been  extensively  studied  as  models  of  combustion  (cf., 
Kapila  (ref  8)),  chemical  reaction  (cf.,  Fife  (ref  9)),  and  population  dynamics  (cf.,  Hoppensteadt 
(ref  10)).  In  addition,  if 

® 

Eq.  (1)  arises  in  convective  flows  (cf.,  Batchelor  (ref  11)).  Therefore,  even  though  this  method 
does  not  address  every  variety  of  problem  possible,  it  would  still  be  a  useful  numerical  tool  in  a 
large  number  of  scientific  areas. 

Many  technological  situations,  including  some  of  those  mentioned  above,  involve  the 
rapid  formation,  propagation,  and  disintegration  of  small-scale  structures.  Some  examples  are 
shock  waves  in  compressible  flows,  shear  layers  in  laminar  and  turbulent  flows,  phase  boundaries 
in  nonequilibrium  processes,  combustion  fronts,  and  classical  boundary  layers.  It  is  this 
increasing  complexity  of  the  physical  problem  that  scientists  and  engineers  want  to  address,  as 
well  as  the  desire  for  reliable  and  robust  software  tools  to  accurately  and  efficiently  describe 
these  phenomena  that  led  to  the  development  of  adaptive  numerical  techniques,  in  general,  and 
AFEM,  in  particular. 

The  adaptive  strategies  utilized  by  AFEM  are:  (1)  error  estimation  coupled  with  (2)  local 
mesh  refinement  (cf.,  e.g.,  Adjerid  and  Flaherty  (ref  12),  Babuska  and  Dorr  (ref  13),  Babuska, 
Zienkiewicz,  Gago,  and  Oliveira  (ref  14),  Bank  and  Weiser  (ref  15),  Berger  and  Oliger  (ref  16), 
Bieterman  and  Babuska  (refs  17,18),  Moore  and  Flaherty  (ref  19),  Shephard,  Qingxian,  and 
Baehmann  (ref  20),  Strouboulis  and  Oden  (ref  21),  Zienkiewicz  and  Zhu  (ref  22),  and  (3)  mesh 
movement  (cf.,  e.g.,  Adjerid  and  Flaherty  (ref  12),  Amey  and  Flaherty  (ref  23),  Bell  and  Shubin 
(ref  24),  Davis  and  Flaherty  (ref  25),  Dorfi  and  Drury  (ref  26),  Dwyer  (ref  27),  Ewing,  Russell, 
and  Wheeler  (ref  28),  Hyman  (ref  29),  Kansa,  Morgan,  and  Morris  (ref  30),  Miller  and  Miller 
(refs  31,32),  Petzold  (ref  33),  Rai  and  Anderson  (ref  34),  Russell  and  Ren  (ref  35),  Saltzman  and 
Brackbill  (ref  36),  Smooke  and  Koszykowski  (ref  37),  Thompson  (ref  38),  Verwer,  Blom, 
Furzeland,  and  Zegeling  (ref  39),  and  White  (ref  40)).  Detailed  descriptions  of  how  AFEM 
implements  these  adaptive  strategies  are  found  in  separate  technical  reports  by  Coyle  and 
Flaherty  entitled:  Adaptive  Finite  Element  Method  II:  Error  Estimation  (ref  41);  Adaptive  Finite 
Element  Method  III:  Mesh  Refinement  (ref  42);  and  Adaptive  Finite  Element  Method  IV:  Mesh 
Movement  (ref  43). 
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This  report  is  an  introduction  to  AFEM  as  well  as  an  illustration  of  AFEM’s  usefulness 
as  a  computational  tool.  In  particular,  the  purpose  of  this  report  is  to  formally  present  the 
equations  that  AFEM  attempts  to  solve,  describe  the  solution  algorithm  employed  by  AFEM, 
and  present  computational  examples.  Specifically,  in  the  Problem  Statement  section,  the 
problem  that  AFEM  was  developed  to  solve  is  formally  stated  and  a  weak  form  of  the  problem  is 
derived.  The  Spatial  Discretization  and  Temporal  Discretization  sections  illustrate  how  AFEM 
numerically  discretizes  the  weak  form  of  the  problem  developed  in  the  second  section  on  a 
nonuniform,  moving  mesh.  The  Spatial  Discretization  section  delineates  the  spatial 
discretization,  demonstrating  how  AFEM  utilizes  both  piecewise  linear  and  quadratic  finite 
elements.  The  Temporal  Discretization  section  recounts  the  temporal  discretization,  depicting 
how  both  the  backward  Euler  method  and  the  trapezoidal  rule  are  applied  to  the  equations 
derived  in  the  third  section.  Then  a  brief  description  of  AFEM’s  adaptive  procedures  is  given 
and  the  overall  solution  algorithm  is  presented.  The  Example  section  presents  computational 
results  for  four  test  problems  chosen  in  order  to  illustrate  Afferent  aspects  of  the  code.  Finally, 
the  last  section  contains  a  summary  and  discussion  of  this  report. 


PROBLEM  STATEMENT 

Consider  the  problem  of  finding  a  numerical  solution  of  an  M-dimensional  system  of 
partial  differential  equations  of  the  form 

subject  to  the  initial  conditions 

u(x,0)  =  H®(x)  ,  a^x^b 

and  linear  separated  boundary  conditions 

A\t)uia,t)  +  =  gXi) 

A\t)uQ>jt)  +  =  g\i) 

The  variables  x  and  t  represent  spatial  and  temporal  coordinates  and  denote  partial 
differentiation  when  they  are  used  as  subscripts; «,/  u°,  and  are  M-vectors;  and  D,  A\  B,  A\ 
and  W  are  MxAf  matrices. 

The  problem  is  assumed  to  be  well-posed  and  parabolic;  thus,  e.g.,  D(x,t,u)  is  positive 
definite.  The  rows  of  B  and  ff  are  restricted  to  be  either  entirely  zero  or  a  row  of  the  Af  xM 
identity  matrix.  When  the  f**  row  of  or  jff’  is  identically  zero,  then  Al  or  .4;-  cannot  be  zero, 
respectively,  and  the  boundary  condition  is  a  Dirichlet  (essential)  condition.  Otherwise,  the 
boundary  condition  is  a  Neumaim  or  Robbins  (natural)  condition. 

A  weak  form  of  the  problem  is  constructed  by  multiplying  Eq.  (3a)  by  a  test  function 
v(x,t)  G  Ho ,  integrating  the  result  with  respect  to  x  from  a  to  b,  and  integrating  the  diffusive 
term  by  parts  to  obtain 


a<x<b  ,  r>0,  (3a) 

(3b) 

,  r>0  .  (^‘^^ 
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(v,Mj)  +  (v/l  +  A(y,u)  =  ,  t>0  ,  for  all  v  £  . 


(4a) 


The  inner  product  (v,u)  and  strain  energy  y4(v,«)  are  defined  as 

(v,if)  =  f^v^udx,  A(v,u)  =  f^’v^Du^  dx  .  (4b, c) 

J  a  ^  c 


Functions  v  belonging  to  are  required  to  have  finite  values  of  (v,v)  and  (v^,vj.  Functions  in 
Hi  are  in  H  and  vanish  at  jc  =  a  and/or  b  if  an  essential  boundary  condition  is  applied  there. 
Any  weak  solution  h  G  of  Eq.  (4a)  must  also  satisfy  any  essential  boundary  conditions  at  j:  = 
a 


or  at  j;  =  fe 


ii.(a,f) 


M 


M 

7=1 


ufbf) 


gft) 


Y,  Alit)ufb,t) 


j*i 

7=1 


(4d) 


(4e) 


when  the  i**  row  of  JB'  and/or  W  is  zero,  respectively.  Natural  boundary  conditions  replace  the  f 
component  otu^zXx  —  a  ox  b  in  Eq.  (4a)  when  prescribed. 


Initial  conditions  for  Eq.  (4a)  are  obtained  by  projection,  i.e., 

(v,u)  =  (v,tt®)  ,  t  =  0  ,for  all  V  E  Hq  . 


(4f) 


A  discrete  version  of  the  weak  system  Eq.  (4)  is  constructed  by  using  finite  element- 
Galerkin  procedures  in  space  and  finite  difference  techniques  in  time  on  a  fully  adaptive  mesh 
(one  that  is  both  refined  and  moved  as  time  progresses). 


SPATIAL  DISCRETIZATION 

To  discretize  Eq.  (4a)  in  space,  introduce  a  time-dependent  partition 

n/f)  =  {a  =  XQ<xft)<X2it)<...<Xff  =  b)  < 

of  (a,b)  into  N  subintervals  xft)),  i=2,^...,Nand  approximate  «  and  v  by  piecewise 

polynomial  functions  V  and  V,  respectively,  with  respect  to  this  partition.  Thus,  the  spatially- 
discrete  form  of  Eq.  (4a)  consists  of  finding  17  G  5^  C  such  that 
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(6a) 


(y,U)  +  +  A{V,U)  =  V'^DUj’,  ,  t>0  , 


for  all  V  e  Sq  <=■  Hq  , 


iy,U)  =  (V,II®)  ,  ^=0  ,  for  all  V  e  c  Hq 


The  spaces  5^  and  5^  will  consist  of  either  piecewise  linear  or  piecewise  quadratic  polynomial 
functions.  The  spaces  of  piecewise  linear  polynomials  are  denoted  and  and  a  basis  is 
easily  constructed  in  terms  of  the  familiar  "hat"  functions 


x-x._ff) 
xft)  -X^_ft) 
4>iix,t)  =  x.,^(f)-x 


,  ^  ^  xff) 


— —  ,  xfi)  ^x  ^ 
0  ,  otherwise 


The  piecewise  linear  finite  element  solution  Vi^  is  written  in  the  form 

N 

u^ipcf)  = 

1=0 


and  determined  by  solving  the  ordinary  differential  system 


+  AiV„U,)  =  vlDU.^l  , 
r>0  ,  for  all  Vj  6  Sq’^  , 


(.V,,U,)  =  (VpU®)  ,  t=0,for  all  V,  e  So‘ 


where  the  piecewise  linear  test  functions  Vj  E  have  a  form  similar  to  Eq.  (8). 

Piecewise  quadratic  approximations  f/2  G  are  constructed  by  adding  a  "hierarchical 
correction  E2(x,t)  to  f/j,  i.e., 

V^{x,t)  =  Uf,x,t)  +  E^ixf)  ,  (10a) 


where 
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(10b) 


N 

E2(x,t)  = 

i=l 


The  basis  rpi.v,(x,t),  i=l,2,...,N,  for  the  quadratic  correction  has  the  form 

2[JC-x._,(0][x-x,(0] 


,  ^  <  xfi) 


0  ,  otherwise 


Piecewise  quadratic  solutions  are  determined  by  solving 


{y^,U^)  +  +  A{V^,U^)  =  vlDU^^  \l , 


t>0  ,for  all  Vj  e  , 


iy^,U^)  =  (y^,u^  ,  t=0  ,for  all  E  , 


(11) 


(12a) 


(12b) 


where  Kj  has  a  form  similar  to  Eq.  (10). 


TEMPORAL  DISCRETIZATION 

Discretization  in  time  is  performed  by  integrating,  for  example,  Eq.  (6a)  over  the  time 
step  from  to  r  to  obtain 


for  all  V  E  Sq  , 


(13a) 


(V,U)  =  (y,H°)  ,  t=0  ,/or  all  V  e 


(13b) 


Rewriting  the  left-hand  side  of  Eq.  (13)  to  emphasize  the  contribution  over  each 
subinterval  gives 
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(14a) 


+  y^/  +  vIdU^  dxdt  =  V^DUjl  dt  , 

for  all  V  €  So  , 


iV,U)  =  (y,tt®)  ,  t=0  ,for  all  V  e  . 


(14b) 


The  integration  in  Eq.  (14)  will  be  over  an  irregular  region  due  to  the  mesh  motion  with  the  test 
function  V  having  an  undesirable  time  dependency. 


In  order  to  overcome  this  difficulty,  introduce  a  linear  transformation 
+  V2Axr'(l+0  +V2(Ax;-Aj:r')T(l  + 


(15a) 


t  =  +  At"T 


where 


Ax,"  =  X,"  -  x,:i  ,  At"  =  t"  - 


tn-\ 


and 


X,"  =  x,(f»)  ,  i  =  0,l,...,Ar 


(15b) 

(15c, d) 

(15e) 


This  transformation  maps  the  rectangle  {(§,t)|-1  <^<1,  0<t<1}  onto  the  trapezoidal 
space-time  element  whose  vertices  are 

and  Oc'f")  ■ 


The  basis  elements  4>i.i(x,t)  and  4>fx,t),  the  only  nonzero  ones  on 
{(x,t)\Ki.ft)  <  X  <  xft)  ,  f  <  t  <  f},  are  transformed  to  functions  with  no  t  dependency;  thus, 
<l>i.fx,t)  and  (f>i(x,t)  become,  respectively. 


=  V2(l-0  ,  and 


(16a) 


=  V2(l+0  ,  -1  :£  5  ^  1  . 


(16b) 


Define 


7 


and  write  Eq.  (14)  as 


N 


E 


V^DUX  dt ,  for  all  V  e  Sq  , 


(18a) 


(V,V)  =  (V;k®)  ,  t-0  ,for  all  V  e  . 

by  substituting  Eq.  (17)  into  Eq.  (14).  Equation  (18a)  can  also  be  rewntten  as 

Y,  F.  =  J*  ,/or  all  V  e  Sq  , 

i=l  ® 


(18b) 


(18c) 


(V,U)  =  (y,«°)  ,  f=0  ,for  all  V  e  , 


(18d) 


by  performing  the  change  of  variables  from  t  to  t  (cf.,  Eq.  (15b)). 

Transforming  Eq.  (17)  from  the  (*:-t)-plane  to  the  (^,T)-plane  (cf.,  Eq.  (15))  yields 


(19a) 


which  reduces  to 


y\ux^\  -  y\ux\f  +  y^fx^f  +  ylou^ 


dldx 


(19b) 


where 


X  = 


(19c) 


Equation  (19b)  can  be  simplified  further  by  noting  that 


(20a) 


=  Af" 

and  integrating  by  parts  to  obtain 

F.  =  G.(l)  -  G.(0)  +  //t)  rfT 

where 

Gf,z)  -  /‘^  d5  , 

If.z)  =  f  *  V%  +  vfWji]  . 

•  ”1 


Substituting  Eq.  (20b)  into  Eq.  (18c)  then  yields 

N 

'  Ig.(1)  -  G.(0)  +  At" 

lo 


£  [g,(1)  -  G^O)  t  Aff  ‘  I,(r)dz 

i=i  L 


=  Af 


’I'  V^OV,t  dz  , 


N 


for  all  V  €  So  , 


(y,U)  =  (y,«®)  ,t=0,forallV€  So  . 


(20b) 

(20c) 

(20d) 

(21a) 

(21b) 


All  that  remains  is  to  approximate  the  time  integrals  in  Eq.  (21)  using  a  quadrature  rule. 
This  is  done  by  using  a  weighted  two-step  method,  which  for  Eq.  (21)  has  the  form 


I;  [G.(l)  +  G.(0)  +  Ar»0/.(1)  +  Af"(l-^/,(0)]  =  At"dV^DuXl._, 
1=1  ^ 

+  Ar"(l-e)V^DGj*l,^,/oraH  V€5(f ,  0G  [0,1] , 


(22a) 


iV,U)  =  (y,u°)  ,  ^=0  ,for  all  V  e  Sq  . 


(22b) 


Only  two  choices  of  0  are  utilized:  either  0  =  1,  which  yields  the  backward  Euler  method,  or  0  = 
Vz,  which  yields  the  trapezoidal  rule. 
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ADAPTIVE  AND  SOLUTION  ALGORITHMS 


The  principle  adaptive  feature  of  AFEM  is  the  incorporation  of  a  mesh  moving  or  an  r- 
refinement  scheme  (r  for  redistribution).  The  basic  idea  of  mesh  movement  is  to  initially  place  a 
fixed  number  of  mesh  points  in  optimal  locations  and  then  move  them  as  the  problem  is  solved 
so  that  their  locations  remain  optimal. 

AFEM’s  mesh  moving  algorithm  is  based  on  equidistribution,  made  popular  by  De  Boor 
(ref  44)  for  problems  in  functional  approximation.  The  goal  of  equidistribution,  at  least  in  the 
static  case  of  functional  analysis,  is  to  determine  a  mesh  so  that  a  particular  quantity  is  uniform 
over  each  subinterval.  In  this  case,  it  has  been  shown  (cf.,  De  Boor  (ref  44))  that  the  task  of 
selecting  a  mesh  to  minimize  the  maximum  discretization  error  is  asymptotically  equivalent  (for  a 
dense  mesh)  to  equidistributing  the  local  discretization  error. 

In  order  to  facilitate  the  use  of  a  static  equidistribution  law  in  a  method  for  solving  time- 
dependent  problems,  a  few  adjustments  are  made.  First,  the  equidistribution  law  is  extended  to 
include  a  time  dependency.  Next,  the  resultant  equation  is  approximately  normalized.  Finally, 
the  normalized  equation  is  differentiated  with  respect  to  time  (cf.,  Coyle  and  Flaherty  (ref  43)). 

The  result  of  the  above  process  is  a  stable  ordinary  differential  equation  which  represents 
a  dynamic  equidistribution  law.  This  dynamic  law  no  longer  demands  that  some  quantity  be 
uniform  over  each  subinterval,  however.  Instead,  a  time-dependent  mesh  is  determined  such  that 
the  same  quantity  as  in  the  static  case  now  remains  constant  for  all  time  over  each  subinterval. 

Mesh  movement  was  the  first  adaptive  technique  developed  for  time-dependent  partial 
differential  equations  because  it  does  not  demand  as  much  computer  resources  as  some  other 
adaptive  techniques.  No  extra  calculations  to  determine  error  approximations  need  be  made  nor 
extra  storage  need  be  created  to  accommodate  additional  mesh  points  and  the  data  structures  are 
simple. 

Mesh  movement  does  have  its  drawbacks,  however.  If  sufficient  care  is  not  exercised, 
mesh  trajectories  can  leave  the  domain  of  definition,  cross  each  other,  or  oscillate  wildly  from 
time  step  to  time  step  (cf.,  Coyle,  Flaherty,  and  Ludwig  (ref  45)). 

Another  drawback  to  mesh  movement  is  that  it  does  not  guarantee  that  the  approximate 
solution  satisfies  a  prescribed  accuracy  tolerance  for  any  discretization  level.  As  a  result,  most 
researchers  have  either  abandoned  moving  in  favor  of  other  adaptive  techniques  that  guarantee 
such  accuracy  or  incorporated  such  methods  into  their  moving  codes,  as  was  done  with  AFEM. 

The  additional  adaptive  technique  incorporated  into  AFEM  is  a  local  or  h-refinement 
method  coupled  with  a  posteriori  error  estimation.  A  local  refinement  method  is  one  that 
determines  subregions  of  the  domain  where  the  solution  has  not  been  adequately  approximated 
and  solves  the  problem  again  using  a  finer  discretization  level  in  these  subregions.  This 
determination  is  usually  made  by  estimating  the  error  of  the  numerical  approximation  throughout 
the  entire  domain.  Local  refinement  methods  for  time-dependent  partial  differential  equations 
do  not  only  refine  the  discretization,  but  also  coarsen  it  in  regions  where  a  fine  discretization  is 
unnecessary.  In  this  way,  the  methods  attempt  to  lessen  the  burden  associated  with  solving  the 
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problem  more  than  once  in  some  region. 

AFEM’s  refinement  strategy  estimates  the  total  discretization  error  as  well  as  its  spatial 
and  temporal  components.  Then  the  type  of  refinement  that  occurs  is  dependent  upon  which 
error  component  is  dominant.  If  the  temporal  error  dominates,  then  a  new  global  time  step  is 
calculated  and  the  problem  is  solved  again  using  this  time  step.  If  the  spatial  component  of  the 
error  is  the  dominant  one,  then  more  mesh  points  are  added  (and/or  deleted)  in  regions  where 
the  estimate  is  high  (and/or  low)  and  the  problem  is  solved  again  with  the  same  time  step  or  a 
larger  one  if  the  temporal  error  estimate  is  sufficiently  small.  If  no  component  dominates,  and 
yet  the  total  error  estimate  is  too  large,  then  refinement  occurs  in  both  space  and  time  (cf.,  Coyle 
and  Flaherty  (ref  42)). 

A  posteriori  error  estimates  are  obtained  by  effectively  solving  the  problem  more  than 
once,  using  higher  order  methods,  and  then  finding  the  difference  between  the  higher  order  and 
the  lower  order  solutions.  Hierarchial  formulation  and  nodal  superconvergence  (cf.,  Adjerid  and 
Flaherty  (ref  12))  are  used  to  reduce  the  computational  complexity  of  the  error  estimation 
procedure  (cf.,  Coyle  and  Flaherty  (ref  41)). 

The  solution  algorithm  for  solving  Eq.  (1)  consists  of  dividing  the  temporal  region  into 
strips  0  =  ^  <  <  ....  Suppose  that  a  solution  has  been  successfully  determined  up  to  some 

time  level,  The  goal,  then,  is  to  compute  a  solution  with  a  given  user  tolerance,  TOL,  in  the 
Hj  norm  at  the  next  time  level  f. 

The  first  step  is  to  advance  the  mesh  using  the  dynamic  equidistribution  law  to  a  time 

level 

fft  ^  ^n-l  +  ^fn-i  (23) 

where  is  the  last  time  step  used  to  successfully  determine  the  solution  at  time  level,  r  ^ 

Next,  Eq.  (1)  is  discretized  in  space  using  Galerkin’s  method  with  piecewise  linear  finite  element 
approximations  (cf.,  Eq.  (9)).  Temporal  discretization  is  performed  by  the  Backward  Euler 
method  (cf.,  Eq.  (22)  with  0  =  1).  A  second  solution  is  calculated  using  trapezoidal  rule 
integration  (cf.,  Eq.  (22)  with  9  =  Vi)  and  the  difference  between  the  two  solutions  is  used  to 
furnish  the  temporal  error  estimate.  A  third  solution  is  obtained  using  quadratic  finite  elements 
(cf.,  Eq.  (12))  and  the  trapezoidal  rule  (cf.,  Eq.  (22)  with  9  =  Vz).  The  difference  between  this 
third  solution  and  the  original  piecewise  linear  Backward  Euler  solution  is  used  to  furnish  an 
estimate  to  the  total  error.  The  difference  between  this  third  solution  and  the  piecewise  linear 
trapezoidal  rule  solution  provides  the  spatial  error  estimate.  The  total  error  estimate  is  now 
compared  to  TOL.  If  the  error  estimate  exceeds  TOL,  spatial  and/or  temporal  h-refinement  is 
performed  until  the  total  error  estimate  at  the  resultant  time  level  is  less  than  TOL. 
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EXAMPLES 


In  this  section  we  examine  the  performance  of  the  adaptive  finite  element  method 
described  in  the  previous  sections  on  four  examples.  In  each  example,  we  try  to  concentrate  on  a 
particular  facet  of  the  method  as  well  as  illustrate  the  method’s  dependence  on  some  different 
input  parameters  available  to  the  user. 

Example  1 

Consider  a  problem  for  Burgers’  equation: 


“t  ^  (24a) 

u{x,0)  =  sin  Tijc  ,  0  ^  jc  ^  1  ,  (24b) 

u(0,t)  =  a(l,0  =  0  ,  f  >  0  .  (24c, d) 

where  e  =  5x  10'^.  The  solution  of  this  problem  is  a  pulse  that  steepens  as  it  travels  to  the  right 
until  it  forms  a  shock  layer  at  =  1.  After  a  time  of  0(e'*),  the  pulse  dissipates  and  the  solution 
decays  to  zero  (cf.,  Figure  la  at  the  end  of  Example  1). 

This  problem  was  solved  for  a  value  of  TOL  equal  to  0.8  but  with  different  initial  values 
of  N  and  At.  In  Table  1,  a  summary  of  the  results  is  presented  when  N  and  At  are  initially  100 
and  0.1,  respectively.  Table  2  presents  a  summary  for  an  initial  N  of  20  and  an  initial  At  of  0.02. 
Table  3  summarizes  results  when  N  and  At  are  initially  50  and  0.1.  Mesh  trajectories 
corresponding  to  the  data  presented  in  Tables  1,  2,  and  3  are  shown  in  Figures  lb,  2,  and  3, 
respectively  (at  the  end  of  Example  1).  The  horizontal  lines  in  these  figures  are  used  to 
distinguish  the  regions  where  the  values  of  At  are  different. 

These  tables  present  the  relevant  refinement  data  for  the  significant  time  levels  of  the 
solution  process.  By  a  "significant  time  level,"  we  mean  a  time  level  at  which  a  refinement  has 
occurred  and  one  immediately  preceding  a  refinement.  The  tables  are  divided  into  four  columns 
labelled  "Significant  Times,"  W  (with  the  initial  N  in  parentheses),  "At"  (with  the  initial  At  in 
parentheses),  and  "Error  Estimate,"  respectively.  The  significant  time  levels  are  found  in  the  first 
column.  Reading  across  the  tables,  one  finds  the  values  of  N  and  At  necessary  to  complete  the 
solution  step  to  within  the  given  tolerance.  A  dash  is  used  to  indicate  no  change  in  AT  or  Ar  from 
the  preceding  row.  The  last  column  lists  the  value  of  the  total  error  estimate  at  the  successful 
completion  of  the  time  step.  The  last  row  of  each  table  contains  results  for  the  final  time  of  0.8. 
In  this  manner,  the  tables  outline  the  refinement  process  giving  some  indication  as  to  why 
refinement  occurred  as  well  as  demonstrating  the  success  of  the  process. 
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Table  1 


Numerical  Parameters  for  Solving  Eqs.  (24)  with  TOL  =  0.8  and  N  and  ^ 
Initially  100  and  0.1,  Respectively 


Significant 

N 

At 

Error 

Times 

(100) 

(0.1) 

Estimate 

0.1000 

- 

- 

0.2272 

0.1619 

- 

0.0619 

0.2161 

0.2238 

- 

- 

0.5470 

0.2713 

164 

0.0475 

0.3692 

0.3664 

- 

- 

0.4848 

0.3970 

- 

0.0362 

0.3693 

0.8000 

- 

- 

0.1606 

Table  2  Numerical  Parameters  for  Solving  Eqs.  (24)  with  TOL  =  0.8  and  N  and  At 
Initially  20  and  0.02,  Respectively 


Significant 

N 

At 

Error 

Times 

(20) 

(0.02) 

Estimate 

0.24 

- 

- 

0.5094 

0.26 

40 

- 

0.2832 

0.30 

- 

- 

0.5806 

0.32 

69 

- 

0.1973 

0.80 

- 

- 

0.2901 

Table  3 


Numerical  Parameters  for  Solving  Eqs.  (24)  with  TOL  =  0.8  and  N  and  At 
Initially  50  and  0.1,  Respectively 


Significant 

N 

At 

Error 

Times 

(50) 

(0.1) 

Estimate 

0.2000 

- 

- 

0.3127 

0.3000 

75 

- 

0.2292 

0.4000 

- 

- 

0.3613 

0.4572 

- 

0.0572 

0.2522 

0.8000 

- 

- 

0.1389 

In  Adaptive  Finite  Element  Method  III  (ref  42),  Eqs.  (24)  are  solved  for  a  value  of  TOL 
=  0.8  and  the  same  initial  values  of  N  and  At  as  above.  Only  h-refinement  was  employed  in 
Reference  42,  however,  since  the  purpose  of  that  report  is  to  introduce  and  discuss  our 
refinement  procedure.  In  Tables  4,  5,  and  6,  the  results  for  solving  Eqs.  (24)  using  only  our 
refinement  scheme  is  presented.  This  is  done  in  order  to  study  the  effects  of  our  mesh  moving 
scheme  when  coupled  with  refinement. 


Table  4  Numerical  Parameters  for  Solving  Eqs.  (24)  with  TOL  =  0.8  and  N  and  At 
Initially  100  and  0.1,  Respectively 


Significant 

Times 

N 

(100) 

At 

(0.1) 

Error 

Estimate 

0.1000 

- 

- 

0.3221 

0.1566 

- 

0.0566 

0.2218 

0.2131 

- 

- 

0.5436 

0.2549 

165 

0.0417 

0.3455 

0.8000 

- 

- 

0.2728 

Table  5  Numerical  Parameters  for  Solving  Eqs.  (24)  with  TOL  =  0.8  and  N  and  At 

Initially  20  and  0.02,  Respectively 


Significant 

N 

At 

Error 

Times 

(20) 

(0.02) 

Estimate 

0.24 

- 

- 

0.6741 

0.26 

46 

- 

0.3710 

0.28 

- 

- 

0.4330 

0.30 

66 

- 

0.3072 

0.38 

- 

- 

0.6210 

0.40 

109 

- 

0.2488 

0.80 

-  . 

- 

0.1595 

Table  6  Numerical  Parameters  for  Solving  Eqs.  (24)  with  TOL  =  0.8  and  N  and  At 
Initially  50  and  0.1,  Respectively 


Significant 

N 

At 

Error 

Times 

(50) 

(0.1) 

Estimate 

0.1000 

- 

- 

0.3263 

0.1573 

62 

0.0573 

0.2227 

0.2147 

- 

- 

0.4766 

0.2566 

81 

0.0419 

0.3913 

0.3405 

- 

- 

0.4662 

0.3824 

118 

- 

0.3448 

0.8000 

- . 

- 

- 

0.1644 

Comparison  of  the  above  tables  is  very  encouraging.  Consider  for  example,  the  case 
where  N  and  At  are  initially  20  and  0.02,  respectively  (cf..  Tables  2  and  5).  In  both  instances,  the 
initial  time  step  was  small  enough  so  that  no  temporal  refinement  was  necessary,  and  the  code 
recognized  this  fact  in  each  instance.  However,  when  moving  was  coupled  with  refinement,  only 
two  refinements  occurred  and  only  69  elements  were  ultimately  needed  to  solve  the  problem 
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within  a  tolerance  of  0.8  (cf.,  Table  2).  For  h-refinement,  three  refinements  were  necessary  and 
109  elements  ultimately  were  used  (cf.,  Table  5).  This  is  a  significant  reduction  in  the  necessary 
level  of  discretization,  as  we  hoped  would  be  accomplished  by  the  incorporation  of  mesh 
movement. 

Comparing  the  case  where  N  and  At  are  initially  50  and  0.1,  respectively,  is  similar  (cf.. 
Tables  3  and  6).  When  moving  was  coupled  with  refinement,  only  two  refinements  occurred, 
only  75  elements  were  ultimately  needed,  and  only  a  time  step  of  0.0572  was  finally  necessary  to 
solve  the  problem  within  a  tolerance  of  0.8  (cf..  Table  3).  For  h-refinement,  three  refinements 
were  necessary  (two  of  them  in  both  space  and  time),  118  elements  ultimately  had  to  be  used, 
and  the  final  time  step  employed  was  0.0419  (cf..  Table  6).  Once  again,  a  significant  reduction  in 
the  necessary  level  of  discretization  was  accomplished  by  the  incorporation  of  mesh  movement. 

It  is  only  in  the  case  where  N  and  At  are  initially  100  and  0.1,  respectively,  that  mesh 
movement  does  not  improve  the  computation  (cf..  Table  1  and  Table  4).  When  moving  was 
coupled  with  refinement,  three  refinements  were  necessary  (one  of  them  in  both  space  and  time), 
164  elements  ultimately  had  to  be  used,  and  the  final  time  step  employed  was  0.0362  (cf.. 

Table  3).  In  the  h-refinement  case,  two  refinements  were  necessary  (one  of  them  in  both  space 
and  time),  165  elements  ultimately  had  to  be  used,  and  the  final  time  step  employed  was  0.0417 
(cf..  Table  6).  No  significant  reduction  in  the  necessary  level  of  discretization  was  accomplished 
by  the  incorporation  of  mesh  movement.  As  a  matter  of  fact,  it  could  be  argued  that  mesh 
movement  performed  worse  since  an  extra  temporal  refinement  step  was  necessary,  while 
employing  only  one  less  element.  This  is  not  discouraging  since,  clearly,  too  many  elements  were 
used  in  the  initial  mesh;  hence,  it  is  difficult  for  a  moving  or  redistribution  scheme  to  make  any 
dramatic  improvements. 
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for  Example  1  corresponding  to  Table  1. 
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Example  2 

Consider  the  following  partial  differential  equation 

+  AU,u,u^  =  l,f>0.  (25a) 

2rf 

The  initial  and  boundary  conditions  and  the  function /are  chosen  such  that 

u{x,t)  =  tanhrj[(jc-l) +r2t]  (25b) 

is  the  exact  solution.  Function  (25b)  is  a  wave  that  travels  in  the  negative  x  direction  when 
and  are  positive  (cf.,  Figure  4  at  the  end  of  Example  2).  The  values  of  and  Tj  determine  the 
steepness  of  the  wave  and  its  propagation  speed.  We  solved  this  problem  for  values  of  and 
equal  to  10  and  1,  respectively,  in  order  to  study  the  effect  of  a  moving  phenomenon  on  our 
adaptive  algorithms.  In  particular,  two  different  values  of  the  error  tolerance  parameter,  TOL, 
were  selected  to  see  if  the  adaptivity  adjusted  accordingly.  An  initial  value  of  N  equal  to  10  and 
At  equal  to  0.1  were  used  in  both  instanees.  Table  7  summarizes  the  results  for  a  value  of  TOL 
equal  to  0.8,  while  Table  8  presents  results  when  TOL  was  set  to  0.5.  Both  tables  are  similar  to 
the  tables  presented  in  Example  1.  Mesh  trajectories  corresponding  to  the  data  presented  in 
Tables  7  and  8  are  shown  in  Figures  5  and  6,  respectively  (at  the  end  of  Example  2).  The 
horizontal  lines  in  these  figures  are  used  to  distinguish  the  regions  where  the  values  of  At  are 
different. 

Table  7  Numerical  Parameters  for  solving  Eqs.  (25)  with  TOL  =  0.8  and  iV  and  At 
Initially  10  and  0.1,  Respectively 


Significant 

N 

At 

Error 

Times 

(10) 

(0.1) 

Estimate 

0.0545 

18 

0.0545 

0.5007 

0.0893 

27 

0.0347 

0.5197 

0.2628 

- 

- 

0.5933 

0.2975 

48 

- 

0.4456 

1.0000 

- 

- 

0.3205 
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Table  8 


Numerical  Parameters  for  Solving  Eqs.  (25)  with  TOL  =  0.5  and  N  and  At 
Initially  10  and  0.1,  Respectively 


Significant 

Times 


Error 

Estimate 


0.0394 


0.0638 


0.0830 


0.2371 


0.2564 


0.4297 


0.4490 


0.6994 


0.7187 


1.0000 


0.0394 


0.0243 


0.3003 


0.3355 


0.3140 


0.3858 


0.2783 


0.4189 


0.2803 


0.4297 


0.1850 


0.3118 


As  expected,  the  smaller  the  tolerance  the  more  work  the  code  must  perform  in  order  to 
determine  an  approximate  solution  to  within  that  tolerance  (cf.,  Tables  7  and  8).  In  particular, 
notice  how  often  the  time  step  is  adjusted  at  the  beginning  of  the  solution  process  for  both 
tolerance  levels.  The  initial  time  step  of  0.1  is  a  rather  coarse  value  and  the  code  does  recognize 
this  fact  and  reacts  accordingly.  However,  a  better  strategy  would  seem  to  be  a  smaller  initial  At. 

Until  this  point,  all  examples  have  been  performed  on  an  adaptive  mesh  that  initially  was 
uniform.  This  is  not  the  only  option  available  to  a  user  of  this  code.  The  software  will  also 
accept  a  user  determined  mesh  or  the  software  will  determine  an  equidistributed  mesh  (cf., 
Adaptive  and  Solution  Algorithms  section)  using  the  initial  data.  The  attempt  that  the  code 
makes  is  to  equidistribute  the  second  spatial  derivative  of  the  initial  condition  since  the  error  one 
makes  by  a  piecewise  linear  approximation  is  proportional  to  that  derivative. 

The  choice  of  initial  mesh  has  an  effect  on  the  adaptive  properties  of  this  method.  In 
order  to  study  this  effect  we  once  again  solved  Eqs.  (25)  with  the  same  values  of  TOL  and  initial 
values  of  N  and  At.  Instead  of  an  initially  uniform  mesh,  however,  we  allowed  the  code  to 
determine  an  equidistributed  one.  The  results  for  TOL  =  0.8  and  TOL  =  0.5  are  presented  in 
Tables  9  and  10,  respectively.  Mesh  trajectories  corresponding  to  the  data  presented  in  Tables  9 
and  10  are  shown  in  Figures  7  and  8,  respectively  (at  the  end  of  Bcample  2). 
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Table  9 


Numerical  Parameters  for  Solving  Eqs.  (25)  on  an  Initially  Equidistributed  Mesh 
with  TOL  =  0.8  and  N  and  Initially  10  and  0.1,  Respectively 


Significant 

Times 


Error 

Estimate 


0.0541 


0.0860 


0.1121 


0.1381 


0.3465 


0.3725 


1.0000 


0.0541 


0.0320 


0.4867 


0.5449 


0.4995 


0.5145 


0.5449 


0.4480 


0.5087 


Table  10  Numerical  Parameters  for  Solving  Eqs.  (25)  on  an  Initially  Equidistributed  Mesh 
with  TOL  =  0.5  and  N  and  At  Initially  10  and  0.1,  Respectively 


Significant 

N 

Error 

Times 

(10) 

(0.1) 

Estimate 

0.0428 

19 

0.0428 

0.3425 

0.0643 

- 

0.0187 

0.3060 

0.3416 

- 

- 

0.4149 

0.3603 

49 

- 

0.2909 

0.7898 

- 

- 

0.4630 

0.8085 

70 

- 

0.1283 

1.0000 

- 

- 

0.2481 

These  results  are  very  interesting.  First,  compare  the  results  for  TOL  equal  to  0.8  (cf., 
Tables  7  and  9).  In  this  case,  it  is  not  clear  which  initial  mesh  strategy  is  best.  One  case  uses 
slightly  less  elements  in  general  (cf..  Table  9),  while  the  other  uses  slightly  larger  time  steps  (cf.. 
Table  7).  One  case  only  performs  spatial  and  temporal  refinements  independently  (cf..  Table  9), 
while  the  other  performs  two  full  refinements  (cf..  Table  7)  but  performs  fewer  refinement  steps 
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overall.  It  seems  that  for  such  a  crude  tolerance,  an  initially  equidistributed  mesh  gains  a  user 
no  real  advantage  and  perhaps  in  some  ways  hinders  the  process. 

On  the  other  hand,  comparing  the  results  for  TOL  equal  to  0.5  leads  to  a  completely 
different  conclusion  (cf..  Tables  8  and  10).  In  this  case,  starting  the  calculations  with  an 
equidistributed  mesh  is  clearly  advantageous.  With  such  an  initial  mesh,  the  code  uses  far  fewer 
elements  (70  versus  103  ultimately)  and  performs  fewer  refinement  steps  overall  (4  versus  6). 
Even  though  this  strategy  produces  a  smaller  ultimate  time  step  (0.0187  versus  0.0193),  the 
difference  is  not  great.  Also,  this  strategy  does  take  a  larger  initial  time  step  (0.0428  versus 
0.0394)  and  determines  its  ultimate  time  step  more  quickly  (one  less  refinement  step). 
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u  =  tanh  10(x-l+t) 


X 


Figure  4.  Graphs  of  Eq.  (25b)  for  selected  instances  of  time. 
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MESH  TRAJECTORIES 
FOR  UNIFORM  MESH 


Figure  6.  Mesh  trajectories  for  Example  2  corresponding  to  Table  8. 
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MESH  TRAJECTORffiS 
FOR  NONUNIFORM  MESH 


Figure  7.  Mesh  trajectories  for  Example  2  corresponding  to  Table  9. 
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Example  3 


Consider  the  same  partial  differential  equation  as  in  Example  2: 


u„,  0  <  X  <  l,t>0. 

XX  ’  ' 

2rf 


(26a) 


This  time  however,  the  initial  and  boundary  conditions  and  the  function  /  are  chosen  such  that 

u(x,t)  =  tanh  Tj  [(jc  - 1 )  +  rjf]  -  tanh  [x  -  (26b) 

is  the  exact  solution.  In  addition  to  a  wave  that  travels  in  the  negative  x  direction  when  r■^  and  rz 
are  positive,  function  (26b)  is  composed  of  a  wave  that  travels  in  the  positive  x  direction  (cf.. 
Figure  9  at  the  end  of  Example  3).  We  solved  this  problem  for  values  of  and  rj  equal  to  10 
and  1,  respectively,  in  order  to  study  the  effect  of  more  than  one  moving  structure  on  our 
adaptive  algorithms. 

There  is  another  facet  of  the  code  that  we  decided  to  test  with  this  particular  problem. 
This  facet  is  related  to  the  dynamic  equidistribution  law  upon  which  our  moving  algorithm  is 
based  (cf.,  the  previous  section  and  Coyle  and  Flaherty  (ref  43)).  Formally  stated,  the 
equidistribution  problem  is  to  determine  a  dynamic  mesh 

n^f)  =  { a  =  jCq  <  < x^{t)  <  ...  < x^.i(0  <Xff  =  b)  (26c) 

so  that 

f KCO  dc  ‘ij"  dc  ,j  =  lA-JV,  ( a  0 .  (2M) 

J  a  a 

where  w(x,t)  >  0,x  G  [a,b],  t  >  0,h  usually  chosen  to  be  a  function  of  the  solution  of  the 
underlying  partial  differential  equation  (cf.,  (ref  43)).  This  equidistribution  problem  has  a 
nonunique  solution  whenever  w(x,t)  ;=  0.  This  difficulty  is  usually  handled  by  imposing  a  lower 
bound  on  w,  e.g.,  it  is  common  to  replace  w(x,t)  by  w(x,t)  +  tj.  This  positive  parameter  tj  not 
only  guarantees  a  unique  solution  to  the  equidistribution  problem,  but  also  can  exercise  a  control 
on  the  mesh  movement. 

Whenever  tj  is  small  compared  to  w,  then  w  will  continue  to  control  the  mesh  movement. 
However,  if  tj  is  large,  then  w  has  little  effect,  the  constant  tj  controls  the  mesh  movement,  and 
hence,  the  movement  is  restricted  since  the  solution  to  Eq.  (26d)  is  no  movement  when  w  is  a 
constant.  Different  values  of  tj  can  be  input  into  the  code  depending  on  how  much  movement 
the  user  wishes  to  allow. 

The  effect  of  the  parameter  tj  can  be  quantified.  Letx,(t),y/r),  and  Uj  denote  mesh 
trajectories  for  r  >  0  such  that 
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J  a  N  ^  ^ 

r^'’w(c,m  =  -f  fVc-odc, 

J  a  N  ^ 


(26e) 


and 


This  implies  that 

p^‘\w(C,t)^r]]dC  =  P%C,O^C+PTit/C, 

J  ^  «i  /X  J  a 


(26f) 


(26g) 


(26h) 


which  further  implies 

/^V?.t)<ic +/*''%«  =  0. 

•’>/*)  •'“y 


(26i) 


Invoking  the  Mean  Value  Theorem  yields 


Xj(t)  = 


. P^y/f)  - 

w(x(f),0+tl  ^ 


- '■ - Uf, 

W(x(f)^)+Tl  ^ 


(26j) 


where  ;cft)  is  some  function  whose  values  lie  between  x/tj  andy/tj.  Hence,  Xyftj  is  a  linear 
combination  of  the  moving  trajectory  y,ftj  and  the  stationary  trajectory  Uj  with  the  value  of  17 
determining  the  dominant  component. 


Furthermore,  Eq.  (26j)  implies 


max  w(x,r) 


(2«k) 


Choosing  tj  to  be  proportional  to  the  maximum  value  of  w  yields 


(261) 


where 
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(26m) 


Ti  =  r  max  w(x,t). 

aixsbfiO 

Two  different  values  of  the  parameter,  17,  were  selected  to  see  how  the  adaptive 
algorithms  reacted.  An  initial  value  of  N  equal  to  20  and  At  equal  to  0.05  were  used  in  both 
instances  with  TOL  =  0.5.  Table  11  summarizes  the  results  for  a  value  of  17  equal  to  10,  while 
Table  12  presents  results  when  77  was  set  to  2.  Both  tables  are  similar  to  the  tables  presented  in 
Example  2.  Mesh  trajectories  corresponding  to  the  data  presented  in  Tables  11  and  12  are 
shown  in  Figures  10  and  11,  respectively  (at  the  end  of  Example  3).  The  horizontal  lines  in  these 
figures  are  used  to  distinguish  the  regions  where  the  values  of  At  are  different. 


Table  11  Numerical  Parameters  for  Solving  Eqs.  (26)  on  an  Initially  Equidistributed  Mesh 
with  TOL  =  0.5  and  N  and  At  Initially  20  and  0.05,  Respectively,  and  17  =  10. 


Significant 

N 

At 

Error 

Times 

(20) 

(0.05) 

Estimate 

0.0360 

41 

0.0360 

0.2595 

0.0655 

- 

0.0295 

0.2945 

0.0950 

68 

- 

0.2841 

0.1540 

- 

- 

0.3268 

0.1836 

102 

- 

0.2929 

0.3017 

- 

- 

0.3507 

0.3312 

135 

- 

0.2759 

1.0000 

- 

- 

0.2025 

32 


Table  12  Numerical  Parameters  for  Solving  Eqs.  (26)  on  an  Initially  Equidistributed  Mesh 
with  TOL  =  0.5  and  N  and  At  Initially  20  and  0.05,  Respectively,  and  rj  =  2 


Significant 

N 

At 

Error 

Times 

(20) 

(0.05) 

Estimate 

0.0382 

39 

0.0382 

0.2623 

0.0764 

- 

- 

0.3243 

0.1147 

71 

- 

0.2794 

0.1911 

- 

- 

0.3278 

0.2233 

112 

0.0322 

0.2659 

0.3845 

- 

- 

0.3512 

0.4068 

- 

0.0223 

0.2666 

0.4515 

- 

- 

0.3458 

0.4654 

- 

0.0139 

0.2505 

1.0000 

- 

- 

0.1032 

The  interpretation  of  these  results  is  not  as  straightforward  as  the  previous  two  examples. 
At  first,  when  the  movement  is  less  restricted  (cf.,  Table  12)  the  code  is  able  to  take  larger  time 
steps  and  use  fewer  elements.  However,  as  the  solution  process  continues,  the  smaller  number  of 
elements  results  in  a  smaller  time  step. 

When  the  movement  is  more  restricted  (cf..  Table  11),  a  reasonably  refined  time  step  is 
determined  quickly,  but  more  and  more  elements  must  be  used  in  order  to  solve  the  problem 
within  the  specified  tolerance.  The  decision  on  the  restriction  of  movement  would  seem  to 
depend  on  what  is  more  costly  for  the  individual  problem  being  addressed.  More  restrictive 
movement  may  result  in  an  ultimately  larger  time  step  but  with,  perhaps,  a  larger  number  of 
elements  being  used  at  each  of  those  time  levels.  Is  that  really  better  than  an  ultimately  smaller 
time  step  with  far  less  elements  and  with  a  succession  of  intermediate  time  steps  being  used 
before  the  final  temporal  refinement?  The  answer  would  seem  to  depend  on  the  magnitude  of 
the  numbers  for  the  problem.  If  the  number  of  elements  that  must  be  used  when  restricting  the 
movement  is  so  large  that  solving  the  problem  at  any  time  level  is  very  expensive  in  terms  of 
computer  time  and  storage,  then  restricting  the  movement  was  not  the  way  to  proceed.  On  the 
other  hand,  if  unrestricted  movement  causes  time  steps  so  small  that  this  is  the  source  of 
unreasonable  computer  times,  then  the  movement  should  be  more  restricted. 
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u  =  tanh  10(x-l  +t)  -  tanh  10(x  -t) 


X 


Figure  9.  Graphs  of  Eq.  (26b)  for  selected  instances  of  time. 
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MESH  TRAJECTORffiS 
FORri  =  10 


Figure  10.  Mesh  trajectories  for  Example  3  corresponding  to  Table  11. 
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MESH  TRAJECTORffiS 
FORr|  =  2 


Figure  11.  Mesh  trajectories  for  Example  3  corresponding  to  Table  12. 
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Consider  an  elastic  impact  problem  for  circular  cylindrical  rods  as  proposed  by  T.  W. 
Wright  (ref  46)  which  satisfies: 

V  =  Wj,  e  =  p  =  u^,  q  =  u^,  (27a,b,c,d) 

v^  =  S^,p^  =  2iQ^-P).  (27g,h) 

Here,  u  and  w  are  dimensionless  radial  and  axial  displacement  components,  p  and  v  are  radial 
and  axial  velocity  components,  e  is  the  axial  strain,  q  is  the  shear  strain,  and  S,  P,  and  Q  are 
axial,  radial,  and  shear  forces,  respectively.  Equations  (27e,f)  are  compatibility  relations,  and 
Eqs.  (27g,h)  are  the  equations  of  motion.  We  also  need  appropriate  constitutive  laws  that  relate 
5,  P,  and  Q  to  e,  u,  and  q,  and  herein  we  simply  use  the  linear  Hooke’s  laws 

S  =  e  +  -^^u,P  =  -^(ye+u),Q^-^-^q,  (27ij,k) 

1-v  1-v  4(1 -v) 

where  v  is  Poisson’s  ratio.  If  v  =  0,  we  have  a  one-dimensional  theory  and  Eqs.  (27)  reduce  to  a 
simple  wave  equation.  However,  if  v  0,  Eqs.  (27)  give  a  two-dimensional  theory  that  is  valid 

when  the  length  L  of  the  cylindrical  rod  is  large  compared  to  its  unit  radius. 

We  use  V  =  0.3,  homogeneous  initial  conditions  and  the  boundary  conditions 

v(0,f)  =0.8919 ,  Q(0,t)  =0 ,  S(ll,f)  =0,  (?(ll,f)  =0.  (271,m,n,o) 

These  conditions  correspond  to  a  cylinder  of  length  L  =  11  that  is  hit  at  r  =  0  by  a  rigid  wall 
that  is  moving  with  a  velocity  of  0.8919.  We  also  added  viscous  terras  eq^^p  ev^  and  ep,^^  to  the 
right  hand  sides  of  Eqs.  (27d,e,f),  respectively,  and  used  a  value  of  e  =  0.01. 

Strictly  speaking,  this  is  not  the  type  of  problem  that  this  method  was  designed  to  solve. 
Equations  (27)  are  a  system  of  hyperbolic  equations  not  parabolic.  However,  the  method  has 
shown  the  ability  to  successfully  solve  problems  even  when  the  diffusion  matrix,  D(x,t,u),  of  the 
model  equations  (cf.,  Eq.  (3)  is  small  (cf..  Examples  1,  2,  and  3).  This  lead  to  the  addition  of  the 
diffusive  terms  as  mentioned  above  and  the  determination  of  whether  such  diffusion  terms  would 
be  enough  to  enable  the  code  to  solve  such  problems. 

It  turned  out,  however,  that  the  introduction  of  artificial  diffusion  was  not  enough  for  this 
particular  problem.  The  backward  Euler  method,  our  basic  time  integrator  is  a  super  stable 
method  that  tends  to  dissipate  oscillations  even  when  they  are  real.  This  would  lead  our  code  to 
conclude  that  the  temporal  error  was  large  and  hence  cut  the  time  step  to  unacceptable  levels. 


To  avoid  this  dilemma,  we  did  not  propagate  the  backward  Euler  solution  when  solving 
this  problem.  Instead,  we  propagated  the  trapezoidal  rule  solution  and  used  the  backward  Euler 
method  to  furnish  the  local  error  estimates.  This  is  not  quite  correct  because  the  backward 
Euler  solution  is  of  a  lower  order  than  the  trapezoidal  rule  one.  However,  other  researchers 
have  obtained  error  estimates  in  this  manner  and  have  named  the  procedure  "reverse"  or  "local 
extrapolation"  (cf.,  Hairer,  Norsett,  and  Wanner  (ref  47)).  It  seemed  a  reasonable  way  to 
proceed  given  our  framework. 

This  problem  was  proposed  to  us  by  a  colleague  (cf.,  Flaherty,  Coyle,  Ludwig,  and  Davis 
(ref  48)).  He  was  mainly  interested  in  solving  for  the  axial  stress  and  comparing  the  numerical 
results  with  his  theoretical  predictions.  He  had  no  real  interest  in  the  other  stress  or 
displacement  components. 

This  code  gives  the  user  the  ability  to  concentrate  on  some  subset  of  the  components  of 
the  vector  system  of  equations  and  ignore  others.  This  is  done  by  having  the  user  base  adaptivity 
on  a  subset  of  the  equations. 

When  we  solved  this  problem,  the  code  was  instructed  to  base  its  adaptivity  only  on  the 
equation  for  the  axial  stress.  Hence,  when  the  tolerance  parameter,  TOL,  was  set  to  0.5,  it  was 
only  the  axial  stress  that  was  determined  to  within  that  tolerance.  Since  the  axial  stress  was 
coupled  to  the  other  components,  we  wished  to  study  the  effect  of  selecting  only  one  component 
on  which  to  base  our  adaptive  algorithms. 

The  results  of  solving  this  problem  are  very  encouraging.  The  refinement  and  error 
estimation  procedures  detected  the  discontinuity  at  x  =  0,  r  =  0  and  immediately  concentrated 
more  elements  there  as  well  as  cut  the  time  step  (cf.,  Figure  12  at  the  end  of  Example  4). 

Initially,  N  was  50  and  At  was  0.01  and  the  code  changed  this  to  iV  =  75  and  At  =  0.0051. 
The  time  step  was  the  changed  to  0.0017  and  no  refinement  was  necessary  after  that.  Experience 
with  this  problem  has  shown  that  no  additional  refinement  was  necessary  more  as  a  result  of  the 
level  of  the  immediate  refinement  than  of  the  mesh  moving.  This  is  probably  due  to  the 
diffieulty  associated  with  the  initial  discontinuity. 

The  axial  stress  was  more  than  adequately  determined  (cf..  Figure  13  at  the  end  of 
Example  4),  as  hoped,  when  compared  to  theoretical  predictions  and  solutions  computed  with  a 
uniformly  fine  discretization  (N  =  400,  At  =  0.0005)  (cf.,  Fipire  14  at  the  end  of  Example  4).  It 
is  known  that  one-dimensional  elastie  waves  are  nondispersive,  whereas  two-dimensional  waves 
are  dispersive.  The  dispersive  nature  of  the  two-dimensional  solution  is  clearly  visible  in  Figure 
13. 


Other  components  were  not  approximated  as  well.  Occasional  oscillations  would  develop 
(cf..  Figure  15  at  the  end  of  Example  4).  However,  the  qualitative  nature  of  these  solutions 
remained  intact  when  compared  to  solutions  determined  using  a  uniformly  fine  discretization  (cf.. 
Figure  16  at  the  end  of  Example  4).  Therefore,  if  a  user  is  only  interested  in  qualitative 
information  on  some  problem  components,  it  seems  possible  to  save  some  computer  effort  by 
instrueting  the  adaptive  algorithms  to  ignore  or  reduce  the  accuracy  of  those  components. 
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Figure  12.  Mesh  trajectories  for  Example  4. 
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SUMMARY 


This  is  the  first  in  a  series  of  reports  whose  overall  purpose  is  to  report  on  the  progress 
of  a  project  whose  aim  was  to  develop  a  reliable  and  efficient  numerical  method  for  solving 
systems  of  parabolic  partial  differential  equations.  In  particular,  the  goal  was  for  the  resulting 
software  to  be  able  to  compute  an  approximate  solution  within  a  user  specified  error  tolerance 
while  using  a  minimum  of  computer  resources. 

The  basic  approach  taken  to  accomplish  this  goal  was  to  combine  the  adaptive  techniques 
of  mesh  movement  and  local  mesh  refinement.  The  reasons  for  this  choice  were  that  mesh 
movement  could  calculate  the  most  reliable  solution  at  a  fixed  discreti2ation  level,  while  local 
mesh  refinement  could  guarantee  an  accurate  solution  for  any  prescribed  error  tolerance.  Local 
mesh  refinement  schemes  can  be  costly,  due  to  the  necessity  of  recomputing  the  solution; 
however,  proper  mesh  motion  should  permit  as  few  levels  of  refinement  as  possible  (cf., 
;feamples  section). 

The  Problem  Statement  section  formally  stated  the  problem  that  AFEM  attempts  to 
solve.  Then,  in  the  Spatial  Discretization  section,  the  discretization  for  both  piecewise  linear  and 
piecewise  quadratic  elements,  defined  on  a  moving  mesh,  was  presented.  The  discretization  for 
two  time  integration  schemes  was  demonstrated  in  the  Temporal  Discretization  section. 

Although  the  finite  element  method  is  now  fairly  common,  the  concept  of  adaptive  finite 
elements  is  still  relatively  new.  An  auxiliary  purpose  of  this  particular  report  was  to  not  only  give 
a  brief  description  of  the  adaptive  and  solution  algorithms  developed  for  use  in  AFEM,  but  also 
to  provide  an  introduction  to  adaptive  strategies  in  general.  This  was  done  in  the  Adaptive  and 
Solution  Algorithms  section. 

In  the  Examples  section,  examples  were  presented  so  as  to  demonstrate  the  workings  of 
the  entire  adaptive  finite  element  code  and  how  the  different  pieces  interact.  The  results 
indicate  that  we  have  satisfied  our  original  objective  to  a  large  degree.  Solutions  are  computed 
to  within  a  prescribed  error  tolerance  with  a  much  coarser  level  of  discretization  than  would  be 
necessary  without  adaptivity.  The  mesh  movement  does,  indeed,  delay  the  necessity  for 
refinement.  Meanwhile,  refinement  does  not  interfere  with  the  stability  and  smoothness  of  the 
mesh  motion. 

In  the  future,  we  intend  to  extend  the  method  to  higher  spatial  dimensions.  The  question 
of  proper  and  stable  mesh  movement  has  retarded  our  progress  in  this  important  area.  With  this 
obstacle  now  overcome  (cf.,  Coyle  and  Flaherty  (ref  43)),  we  feel  we  can  more  than  adequately 
address  this  topic. 
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